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We consider a problem of nuclear waste contamination. It takes into account the thermal effects. The 

1 <~| temperature and the contaminant's concentration fulfill convection-diffusion-reaction equations. The ve- 
locity and the pressure in the flow satisfy the Darcy equation, with a viscosity depending on both con- 
centration and temperature. The equations are nonlinear and strongly coupled. Using both finite volume 
and nonconforming finite element methods, we introduce a scheme adapted to this problem. We prove 
the stability and convergence of this scheme and give some error estimates. 
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1. Introduction 
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A part of the high-level nuclear waste is now stored in environmentally safe locations. One has to 
consider the eventuality of a leakage through the engineering and geological barriers. It may cause 
the contamination of underground water sources far away from the original repository's location. In the 
\ present paper, we consider such a problem of nuclear waste contamination in the basement. We take into 

account the thermal effects. The evolution in time of the temperature and of the contaminants concen- 
tration is then governed by convection-diffusion-reaction equations. The velocity and the pressure in the 
flow satisfy the Darcy equation, with a viscosity depending on both concentrations and temperature in 
a nonlinear way. The velocity satisfies an incompressibility constraint. We introduce a scheme adapted 
to this problem. We use both finite volume and nonconforming finite element methods. It ensures that 
a maximum principle holds and that the associated linear systems have good-conditioned matrices. We 
prove the stability and convergence of the scheme and give some error estimates. 

Let us briefly point out some previous works. A complete model coupling concentrations and pres- 
sure equations is very seldom studied, since the system is strongly coupled. Instead, each equation 
is considered separately. In the general context of convection-diffusion-reaction equations, numerous 
schemes are available (see \15] or lfl9ll and the references therein). Finite difference schemes are some- 



>< 



times used for the convective term (in 1 20] for instance). But they are not adapted for the complex 
geometry of a reservoir. More recently, finite volume methods were developed and analysed. Let us just 
cite the book |01 or 1 13 ] and the references therein. Finite elements (for the diffusive term) and finite 



volumes (for the convective term) are coupled for instance in 

Hi 

. In convection dominant problems, 
the equations are of degenerate parabolic type. This setting is considered in |@, ITill . The reaction 
terms are specifically studied thr oug h operator splitting methods in Il4ll . Now, in the specific context of 



porous media flow, we mention 11711 who consider only the evolution of the pressure. In III 111 1 211 a more 



complete set of equations is used, and a mixed finite element approximation is developed. We stress that 
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in all theses works, as in most, the thermic effects are not taken into account. 

The present paper is organized as follows. Section [2] is devoted to the derivation of the model. In 
section[3] we introduce the discrete tools used in this paper. It allows us to define the numerical scheme 
of section|5] The analysis of the scheme uses the properties of section@] We then prove the stability and 
convergence of the scheme, in sections [6] and [7] respectively. We conclude with some error estimates in 
section[8] 



2. Model of contamination 

The thickness of the medium is significantly smaller than its length and width. Hence it is reasonable to 
average the medium properties vertically and to describe the far-field repository by a polyhedral domain 
Q. of K 2 with a smooth boundary d£l. It is characterized by a porosity and a permeability tensor K. 
The time interval of interest is [0, T\. We denote by p the pressure, by (a/)^i the concentrations of the 
N r radionuclides involved in the flow and by 9 the temperature. The Darcy velocity is represented by 
u. We assume a miscible and incompressible displacement. Due to the mass and energy conservation 
prin ciples, the flow is governed by the following system satisfied in £2 x [0, T], with i = 1, ..,N r — 1 (see 
El). 

N r — 1 

tyRid,cii + div(a,u) — div(0Z) c Va,-) = s,- — say — A,/?,- (j) a,- + ^ kjj XjRj (j) aj , (2.1) 

0#p^0+div(0u)-div(0D e V0) = -j fl - j (0 - 0*) , (2.2) 

divu = s, u+ Vp=f. (2.3) 

,0) 

In ( 12. II ) the retardation factors Ri > are due to the sorption mechanism. The real A,- -1 > denotes the 
half life time of radionuclide ;. The term —A, /?, a, describes the radioactive decay of the i-th specy. 
Meanwhile, the quantity Y*j^ikijkjRj§aj is created by radioactive filiation. The molecular diffusion 
effects are given by the coefficient D c > 0. The contamination is represented by the source term Sj and 
s = J^lj Sj. In ( 12.21 ) the coefficient c € p > is the relative specific heat of the porous medium. The thermic 
diffusion coefficient is denoted by Dq > 0. The real 0* > is a reference temperature. The constitutive 
relation ( 12.3b is the Darcy law and f is a density of body forces. For a large range of temperatures ji has 
the form 

where fig is a nonlinear function. For instance, in the Koval model for a two-species mixture lfl6ll . we 
have 

Alfl (fl)=M(0)(l + (M 1 / 4 -l) fll )- 4 
where M = ju(0)/jll(l) is the mobility ratio. 

We notice that the equations ( 12. Il l- ( 12731 ) are strongly coupled. Moreover, every concentration equation 
( 12.11 ) involves a different time scale. Therefore, it is difficult to build a numerical scheme that captures 
all the physical phenomena. We have to transform these equations. We first assume that only serial or 
parallel first-order reactions occur, so that kij = y, with y\ = 0. Next, following J3l, we assume that no 
two isotopes have identical decay rates and we set 



d=a u c i = a i + 1 £ [YItZT; ) a J &»i = 2,..,N r . (2.4) 



I A, -A, 



7=1 \l=j 
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Lastly, without losing any mathematical difficulty (see remark l2l2l below). we set Rj = 1 for i=l,..,N r — 
1 and0 = 1, % = 1. We also set s Ci = s, (YYiZ) fr^) Sj and K{c,d) = K/ ^{{ ai )^\d). The 

contamination problem is now modelized by the following parabolic-elliptic system 



d t Ci + div(c,u) —D c Acj = s Cj —scj — A,c ; - , (2.5) 

<9,0+div(0u)-D e 40 = -s -s(d-d*), (2.6) 

divu = s, U+K(c,8)Vp=f, (2.7) 

with i = l,..,N r — 1. These equations are completed with the boundary and initial conditions 

Vc r v = 0, V0v = O, un = 0, (2.8) 

Q| f =o = 4, e\ t=0 = eo. (2.9) 



The pressure p is normalized by f n pdx = 0. Equations d2.51 l are all similar. Thus, for the sake of 
simplicity, we will assume that there is only one. We set N r = 2 and c = c\, Co = c\, s c = s ci , X=X\. 
The results of this paper readily extend to the general case. 

We conclude with some notations and hypothesis. Let D be a bounded open set of with k ^ 1 . 
We denote by the set of functions that are continuous on D together with all their derivatives, 

and have a compact support in D. For p £ {2,°°}, we use the Lebesgue spaces (L P (D), || .||l/>(d)) 
and (L' 5 (D),||.|| L /> (D) ) with U = {LP) 2 . We also use the Sobolev spaces W p ' q (D) for p e and 
q £ In the case D = Q we use the following conventions. We drop the domain dependancy. 

We denote by |.| (resp. ||.||«>) the norms associated to L 2 = L 2 (£2) and L 2 = L 2 (Q) (resp. L°° and 
L°°). We set h\ = {v e L 2 ; J n v(x)dx = 0}. For p G [l,°o[ we define (H p ,\\.\\ p ) and (H p ,\\.\\ p ) with 
H p = W pl an( j hp = (H p ) 2 . Now let (X, |.|) be a Banach space. The functions g : [0, T] -> X such 
that f \\g(t)\\x is continuous (resp. bounded and square integrable) form the set ff(0,T;X) (resp. 

L°°(0,T;X) and L 2 (0,T;X)). The associated norm for the space U°{Q,T;X) (resp. L 2 (<d,T;X)) is 

/ r \l/2 

defined by \\g\\ L ^ { Qj- X ) = sup re[0 , r] ||g(f)IU (resp. \\g\\ 

L 2 (0,T;X) — [Jo 11/(0 III* )■ Finally in all 
computations we use C > as a generic constant. It depends only on the data of the problem. 
We assume the following regularities for the data in ( 12.5I )-( |2~7| > 

KeW 1 '°°((0,l)x(0,°°))> s^^sgeL 2 , fe^(0,r ; L 2 ). (2.10) 

We also assume that k ^ Kj„f with jq„/ > 0. For the initial data, we assume that cq £ H l , 8o & H l , and 
that we have a.e. in Q 

0<c (x)<l, 0_ < o (x) s: 9+ (2.11) 

with 0_ > 0. Finally we assume that we have a.e. in Q 

2s(x)+X^s c {x)^0, 2(0- -9*)s{x)+s e (x) ^0, 2 (0+ - 0*) s(x) +s e (x) ^ 0. (2.12) 

These conditions ensure a maximum principle (proposition ^, ll below). 

REMARK 2.1 We have assumed that first-order reactions occur, so that the coefficients kij in ( 12. Il l 
depend only on i. If fey depends on / and j, one can still uncouple the equations by iterating the transfor- 
mation ( 12.41 . provided that k\j = for j = 2,..,N r — \. This assumption means that the first long-lasting 
isotope disappears and is not created anymore. It is satisfied by many radionuclides. 
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Remark 2.2 We have assumed that the retardation factors Rj are identical. If it is not the case, the 
difficulty and the approach remain the same. Indeed, let us consider the Fourier transform of ( 12. U . For 
a Fourier mode dj(k,t) we obtain 

— (<j)Rj&j) = -<j)(XjRj + lc 2 D c — ikv^aj+yjRj-itytij-i = -Aj/? / -0a 7 +3' I /? 7 -i0fl/-i 

with Xj = Xj + (k 2 D c — iku)/Rj for j = 1, ..,N r — 1. A transform analogous to ( 12.4b uncouples the 
problem. By taking the partial differential equation counterpart, we obtain an equation similar to J2.51 l. 

3. Discrete tools 

3.1 Mesh and discrete spaces 

Let 27~h be a triangular mesh of £2 . The circumscribed circle of a triangle K £ ,% is centered at xk and 
has the diameter Iik- We set h = max^g ^ Hk- We assume that all the interior angles of the triangles of 
the mesh are less than f, so that £ K. The set of the edges of the triangle K £ ,% is S^. The symbol 
n^ i<T denotes the unit normal vector to an edge a £ Sjt and pointing outward K. We denote by (f/, the set 
of the edges of the mesh. We distinguish the subset S'l"' C <§h (resp. S^ ) of the edges located inside £2 
(resp. on dQ). The middle of an edge a £ is x CT and its length is \a\. For each edge a £ let 
and L a be the two triangles having a in common; we set d a = d(x.K a ,XL a ). For all a £ only the 
triangle K a located inside Q is defined and we set d a = d(xK a ,x a ). Then for all <7 £ $h we set T a = 



We assume the following on the mesh (see [13]): there exists C > such that 

\/C7£4, a^CH and |<j|>C7z. 

It implies that there exists C > such that 

Voe^f, t ct = |ct|/^ ct >C. (3.1) 
We define on the mesh the following spaces. The usual space for finite volume schemes is 

Pq = {q £ L 2 ; VA" £ is a constant}. 

For any function qi, £ Pq and any A" £ ,% we set ^ = <27i|a:- We also consider 
/f = {qeL 2 ; V^£ 5/„ q\ K is affine} , 
= £ M > 9a i s continuous over 12} , 
= {q h ePf ; Vcr £ J}"" , g/, is continuous at the middle of a}. 

We have P{ C Z/ 1 . On the other hand P'{ c <£H\ but C //J with 

//j = {q £ L 2 ; VA" £ q\ K £ //'(A")}. 
Thus we define : Hi — » L 2 by setting 

Vq h €H}, VKe.%, V h q h \ K = Vq h \ K (3.2) 
and the associated norm j| . || i /, is given by 

Vq h £Hl k/ I !|i./ I = (b| 2 + |V, I ^| 2 ) !/2 . 
We then have the following Poincare-like inequality for the space P" c HLq (see 12]). 
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Proposition 3.1 There exists C > such that \q h \ < C\V h q h \ for all q h G P" c C)Lq. 
We also define discrete analogues of the norms H l and H~ l for the space Po by setting 

I 1/2 A II II (lh,Wh) 



\9h\\h 



( E T o {qia - qKaf) and ||^ /l ||_ 1)A = sup 



" II II 
V h eP \\Wh\\h 



for any function q/ 2 G Po. Note that for any pi, G Po and g/, G Po, (ph,<lh) ^\\<lh\\-i ,h hh\\h- The following 



Poincare-like inequality holds for the space Pq Lq (see 111310 . 

Proposition 3.2 There exists C > such that \q h \ < C for all # A G Po HLq- 

Finally we set Po = (Pq) 2 , Pf = (Pf ) 2 and use the Raviart-Thomas spaces (0] 

RTjJ = {v/, G Pf ; VKe eS'k, v>,\k- kk.o is a constant} , 

RT = {v,, G RTJ] ; Vtr e 4" r , v h \ K<r -n K<r , a =v,,\ L<y -n K ^ a and v /! -n|^=0}. 

For all \ h G RT , G -5/, and <7 G #x-, we set (v/, • n^^cr = v/,|a" • n^. 
3.2 Projection operators 

We associate with the spaces of section [3~T1 some projection operators. First, we define TIpc : H\ — > P[ 
by setting 

VgGHi, V^GPf, (V(np f9 ),V^) = (V 9 ,V^). (3.3) 

Next, we consider the space Po- Let % = {qv, G L 2 ; q/, is equal a.e. to a continuous function}. We 
define TIp Q :L 2 — ► Po and il/> : % — ► Po by setting 

{^PoP)k = tL / p(x)rfx, (JT/fotfk = , (3.4) 



AT ^ 



for all p G L , q G % and K <E We also set i7p = (iT/> ) . For the space P" c , we define TIp^- : L 
f and n P -c : H l -> P\ 



P" c and Ilpnc : H l -> P," c . For all p G L 2 and q G TV p and JT^ satisfy 



V vffc G P; ,e , (17^/7, i^) = (p, y/ h ) ; Vcr G S h , / (il^?) da = / 9 </<r. (3.5) 

Jo Jo 

For the space RTo, we define iTRT : L 2 — > RTo and nRx : H 1 — > RTo. For all v G L 2 and w G H 1 , 
^RTo v and -TZrtjW satisfy 



V W/I GRT , (n RTo v,w ft ) = (v ) w ft ); Vastf?, / (w-iT RTo w) -n^da = 0. (3.6) 

The operators TIp , Tlpnc (resp. rip , ITrtj) are L 2 (resp. L 2 ) projection operators. They are stable for 

the L 2 (resp. L 2 ) norms. The operators Ilpg, TIpnc and TZrt are interpolation operators. The following 
estimates are classical (Jfl p. 109 and H). 

Proposition 3.3 There exists C > such that for all qeH 1 and v G H 1 

<CA||?||i, |v-n RT „v| <CA||v||i. 
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For all p G H l and q G i/ 2 we have 

\ P - npnc P \ c/2 ||p|| , , \v h ( q - npncq)\ <: ch \\ q \\ 2 . 

For all q G Hf we have 

l^-n^Kc/zll^ll^. 

Finally, using the Sobolev embedding theorem, one checks that 

\np q-n Po q\^.Ch\\q\\ w i,r (3.7) 

for all qeW Lr with r > 1 (see JH). 
3.3 Discrete operators 

Equations ( 12.5b — ( 12.7b use the differential operators gradient, divergence and laplacian. We have to define 
analogous operators in the discrete setting. The discrete gradient operator V/, : P" c — > Pq is the restriction 
to P" c of the operator V/, given by ( 13.2b . The discrete divergence operator div/, : Po — > P" c is defined by 

Vcr G 4"', (div* v„)( X<T ) = . 3|g ' , (v Lct - vjfj • n K , a , 



VciG^r, (div„v ft )(x CT ) = - 



+ IV 

v K a ' n K,a, 



K L , 

for all v/, G Po- It is adjoint to V/, (proposition ^, ll below). The discrete laplacian operator A/, : Pq — > Po 
is the usual one for finite volume schemes (see IU3I1 ). For all q^ G Pq and K E we have 

= t4t E ^ v?At - ^a) • (3-8) 

Let us now consider the convection terms in ( 12.5b and ( 12.6b . We define b : H 1 x H l — > L 2 by 

b(v,*)=div($v) (3.9) 

for all q £ H 1 and v G H 1 . In order to define a discrete counterpart to b we use the classical upwind 
scheme (see IU3I1 ). The discrete operator b/, : RTo x Po ^ Po is such that 

H(Vh,qh)\ K = E 1^1 ((vft-n J r, CT )+^ + (v ft -njc i(T );^z. <J ) (3.10) 

for all v/ ; G RTo, G Po and K G =5/,. We have set a + = max(a,0) and a~ = min(a,0) for all a G KL 
Integrating by parts the convection terms also leads to consider b : L 2 x L 2 x U° — > R defined by 



b(y,p,q) = — p\-Vqdx (3.11) 

for all v G L 2 , p G L 2 and q £ Lf°. The discrete counterpart is b/, : RTo x Po x Po — ► R with 

b h (vh,Ph,<lh) = E ^ x E l CT l (( v /i- n «,<T)ip^ + (v/i-njr, (T )CTPL (J ) (3.12) 

for all v/, G RT , p/, G Po and G Po. 
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4. Properties of the discrete operators 

The properties of the discrete operators are analogous to the ones satisfied by their continuous counter- 
part. The gradient and divergence operators are adjoint. For the operators V/, and div/, we state in 112 111 
the following. 

Proposition 4.1 For all \ h G Po and q h e P" c we have (v^, V/,^/,) = -(<?/,, div^v/,). 
Let us now consider the convection terms. Let q G L°° Hi? 1 , v G L 2 with divv G L? and divv(x) ^ 



a.e. in X2. We obtain b(\,q,q) = f n (q 2 /2) div \dx ^ Oby integration by parts. Forb/, we state in 12 If a 
similar result. 

Proposition 4.2 Let v/ ? G RT with divv/, ^ 0. We have bh(vh,qh,qh) ^ for all q h G Po ■ 

The following stability properties are used to prove the error estimates in section [8] 

Proposition 4.3 There exists C > such that for all p h G Po, qh G Po and v/, G RT with divv/, = 

\H{Vh,Ph,qh)\ ^C\v h \\\p h \\ h \\q h \\ h . (4.1) 
There exists C > such that for all G Po, g/, G Po H Lq, v/, G RTo 

\bh(Vh,Ph,qh)\ ^C(\v h \ ||p/ ? !|/, + |divv/ ? | \\ph\\oo)\\qh\\h- (4.2) 
Proof. For all K G and ae^n we write 

(v/, ^^,0-)^^+ (v/, •n/f, a ) a pL CT = (v/ 7 •n/r. (T ) (J /7A:- \{yh-^K,a)a\ {pL a -pa)- 
Thus (|3.12t reads b/^v/,,/?/,,^/,) = Si +S2 with 

Si = - £ £ I "! I(vft •n K)tT ) tT | (pz, CT , S 2 = P K( 1 K L l ty l( v A- n «,(r)CT. 

Using the Cauchy-Schwarz inequality we write 

E \°\ |va(x ct ) -njf^l {pl -Pk) {qL a ~qK) 

0||v*||-( £ (PL n -^) 2 ) 1/2 ( I fe CT -^J 2 ) 1/2 . 

Since v/, G RT C (Pf ) 2 we have /; 1 1 v/, 1 1 «, ^ C \ \ h | (0] p. 112). Moreover £□) implies "L a e<sj"> (PL a ~ 
PK a ) 2 < C I, ae ejn, To {pl - PK a f =C\\p h \\l and I CTG(f ,, w (g i<7 - qK a ) 2 ^ C \\q h \\\. Thus 

|5i| ^C\\ h \\\p h \\ h \\q h \\ h . (4.3) 
We now consider S<i- We have (v/, • a ) a = for all K E and a G <£a" H ^/f*- Thus we write 

£ |o , |(vA-njt, ff )a= £ l ff K v *' n ^)"= / diWftdx. 

o-G^n*?,;'' 1 ' o-G<sk Ja: 
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It gives the following relation. 

S2= P K % K I divv/,dx = / Pi,qh<tiw h dx. 

Thus if div v/, = then S2 = and estimate ( 14.31 ) gives d4.U . Let us prove d4.21 i. Since £ Pq n Lq we 
can apply proposition [372] Using the Cauchy-Schwarz inequality we get 

\Si\ < ||/»a||~ |divv A | <C||/> A ||oo||^ A || A |divVft|. 



This latter estimate together with ( 14.31 ) gives ( 14.2 
Lastly, we c 

Propositk 
v-n|^ =0 



Lastly, we claim that b/, is a consistent approximation of b [21]. 

Proposition 4.4 Let r > 0. There exists C > such that for all functions q E H 2 and v € H 1+r with 



II^K v >2)-M^RToV,.n/> #)||_ M < cA||?||i||v||i +r . 

Let us now consider the discrete laplacian. We have a coercivity and stability result. 
Proposition 4.5 For all pf, e Po and e Po, we have 

-(A h pi„Ph) = \\ph\\l, \(AhPh,qh)\ < IbftlklkftlU- 

Proof. Definition ( 13.8b implies 

{A h PhAh) = £ L T <y{pL a -pK <r ) = - £ T o{pL a -pK c ){qL a -qK a ). (4.4) 

Setting qi, = pi, gives the first part of the result. Using the Cauchy-Schwarz inequality, we get the second 



We also deduce from ( 14.41 ) the following property. 

Proposition 4.6 For all ph € Pq and q h e fb we have (A h p h ,q h ) — (p h ,A h q h ). 

Lastly, we state that A/, is a consistent approximation of the laplacian. The proof follows the lines of the 



one of proposition 1.14 in 12211 . 

Proposition 4.7 There exists C > such that for all q&H 2 with Vq ■ n|^ = we have 

\\n Po {Aq) - A h (n Po q)\\_ x . h ^Ch\\q\\ 2 . 

5. The finite element-finite volume scheme 

We now introduce the scheme for (l2.5b - d2T9b . The interval [0,7'] is split with a constant time step 
k = T/N. We set [0, T] = U^ofcuWi] with t m = tnk. The time derivatives are approximated using a 
first order Euler scheme. The convection terms are discretized semi-implicitly in time and the other ones 
in an implicit way. We set Sf, = Ilp s, s c h = Ilp s c , s e h = Ilp sg and fj" = f7p f(f m ) for all m £ {0, . . . ,N}. 
Since flp (resp. iTp ) is stable for the L 2 (resp. L 2 ) norm we have 

\s h \^\s\, K\ <|, c |, \s e h \<\s e \, K\ < \f(t m )\ < llfll^LY (5-1) 
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The initial values are = TIp q cq and 9® — Ilp 9o. Then for all n G {0, ... ,N — 1}, the quantities 
c'l + 1 G Pq, 9;'+ 1 G P Q , p" h + 1 G P[ K n L§, u;;+ 1 G RTo are the solutions of the following problem. 

^ -D c A h 4 +1 =4- {s h + X)c n + 1 -b h (u n h ,c» +1 ), (5.2) 

d '^ k e " -D e A h 9;; +1 = -sf 1 -s h (9 n h + l -9,)-b ll (ul 7 9' i ; +l ), (5.3) 

div* ( < +1 v hP ';+ l ) = diva ff 1 - ^ , (5.4) 
u;; +1 =n RTo (^ +1 -< +1 v /7/ ,;;+ 1 ), (5.5) 

with K^ +l = K(c" h +1 , G Po- This term is defined thanks to proposition 16.11 below. Note also that 

the boundary conditions are implicitly included in the definition of the discrete operators ( section [331 ). 
The existence of a unique solution to ( 15.2b and ( 15.3b is classical (see 11311 ). Since K^ +1 ^ fc m ,„ > and 
p1 +l G l\ equation ( 15.41 ) also has a unique solution (see |@]). We have a discrete equivalent for the 
divergence condition (12.7b . 

Proposition 5.1 For all m G {1, . . . ,N} we have divu™ = s h . 

Proof. Let m G {1, . . . ,N} and n = m — 1. We compare the solution of ( 15.4b - d5"31 ) with the solution 
of the following mixed hybrid problem. Let Sq = {jx h : U ff6 ^ — ► M; VtT G S'h, Hi,\a is constant}. Then 
a;; 1 G RT|}, p£ G P and A A m G So are the solution of (see 01) 

Vv„eRTg, «,va) + £ k£ £ |a'| A™, (v/,k -n^) - £ |tf|j«divvak = (f?,ir fc ) , (5.6) 



VjUaG^o, E / M/,(u/"-n)c/c7 = 0, VK€&h, divvtf dx= sdx, (5.7) 

and ^ G is defined by ! a p",lda = A™ for all a G 4. Let ae4 We define </> CT G by setting 
0a(xa) = 1 and <M x a') = for all a' G We set v A = V,,^ eP C RT[J in (gS). We have 

£ Kf £ |a'|A^Va^k-n^= £ KfVa^k £ |a'|A>^ 

and according to the gradient formula 

E |a'|A>*,a'= E / Ph^,cda'= [ V h pjdx. 

&eS K o'e<S K J& Jk 

Thus we get from (15.6b 

(11?, V^) + (KjTV^.V^a) = (C V A a ). (5.8) 
The first term in ( 15.8b is treated as follows. Integrating by parts we get 

(tta n ,V ft <M = -(0 CT) divua n )+ I E / -MuJ'/k-n^cr'. 

Since ( 15.7b implies that u™ £ RTo, we have 



E E ll>cK-n K<& )da'= E |<r|M«a)(a*k-n«a,a-O?k-nxi„a) = 0. 

Ke,%&eS K Ja ' aes'"' 
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Thus (u™, VhQa) = — (0ff)divu™). Then, using (15.71 ). we get 

Furthermore, according to proposition 14.11 we have (k^'V/,/^ 1 , V/,0 ct ) = — (^div/^K"™ V/, /?/")) and 
(f™, V/,0 CT ) = — (0 CT ,diVftf^). Hence we deduce from ( 15.8b that 

V0 ff g /f , (0 CT ,div A «v^) - divftf,; 1 + n^) = o. 

Since {ty a )aeg h is a basis of P" c , we get div/ ! (K - ™V /l /?^) = div/,f™ - Ilpncs™. Thus, by ( 15.4b . there exists 
a real C such that = pf + C. We now compare u™ with u™. Since for all v/, £ RTo we have 

£ < £ |o / |A5(v*k.n i y)= £ | t 7|^(x <T )(u^ CT -n^, CT -<k CT -n JfCT)CT ) = 0, 

/re 55, ct'g^ ctg^' 

it follows from dH that (fi£,v) = (f;; 1 - ^V^.v) for any v G RT . It means that 

u A = n RTo (f/, - k- /; V/,p /? ) = n RTo (t h - k,, w hPh )=u h . 

Thus uj" = u" ; ' satisfies (15.7b and divu™ = ■ 
6. Stability analysis 

We first check that a maximum principle holds. It ensures that the computed concentration and temper- 
ature are physically relevant. 

Proposition 6.1 For any m G {0, . . . ,N} we have ^ < 1 and 0_ 9™ < + . 

PROOF. We prove the result by induction. Since c9 = f7/> co and 0^ = Ilp 9o the result holds for m = 

thanks to ( 12.1 lb and ( 13.41 ). Let us assume that it is true for m =n G {0, ,iV — 1}. Let K € i^,. Equation 

( 15.21 ) implies 

(i+^+U)4 +1 = 4+£4+£A. £ *a(c£ 1 -4 + V*b*(uX > c2 +1 )| jr 

We consider the last term of this relation. Since for any cr G <£fcn <^"' we have 

I'll" n "1+^."+! i ^ ^,"+1 ('.." ^ ^"+1 I / ,,n „ \ + IJl+l „n+\ \ 

{XX h -TiK,a)a c K +{ u h- n K,a)a c L a = {^h ' n K,o)o C K + (-U h ■ n K .o) a (C K -C L<j ). 

We deduce from ( 13.10b 

-b /j K, c ;: +I )| Jf = ^(-4 +1 £ |<T|(uJ.n, )ff ) ff + £ (-uX-n^+^-c^ 1 )). 

Since u£ G RT , (uj • n^)^ = for any ff6 4n 4°*. It implies that E ff6&n<? w |<x| (ujj ■ n K , a ) a = 
T.oe£/c l°1 ( u /l ' n AT,(j)(j- Thus using the divergence formula and proposition 15 . 1 1 we obtain 

1 If 
mFT E l ff IW- n ^V= / divuldx = s K . 

1^1 aeftrW- 11 
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Therefore we get 



{\+2ks K + kX)c"+ i =c" K + ks c K + kD c £ *a(c£ l -c n K +1 ) 



± £ (-u^-n^+^-cf 1 ). (6.1) 



We consider /f, e ^ such that c^ +1 = min^g^ c^ +1 . According to hypothesis ( 12.111 ) and definition 
d3.4t we have 2sR t + A ^ and s^, 0. Thus, using the induction hypothesis, we deduce from ( 16. U 

min c" +I = c" +1 ^ c K i +ks K i J* Q 

We now consider X, € ,% such that cj^~ = max/jg,^ c^. +1 . Using again hypothesis d2.1 II ) we have 
2i^ 5 + A ^ s^- ^ 0. Thus, using the induction hypothesis, we deduce from d6.U 

ma X 4 +1 = C " +1 < fVt^ <; - ^ 1. 

tfe,^ x 1+2^+ArA l+2ks Ks +kX 

A similar work for equation ( 15.3b proves that 0_ ^ min^ G ^ and max^ e « ^ + . Thus the 
induction hypothesis still holds for m = n + 1 . ■ 
We now state the stability of the scheme (I5.2b - d53b . 



Proposition 6.2 For any 1 ^ m < N we have 

m m 

ktHWh+ktimil^C, (6.2) 

n= 1 n= 1 

K| + |V A ^|<C. (6.3) 
PROOF. Let ^ n ^ iV- 1. Multiplying (IBT2l i by 2A:c)^ +1 we get 

«+ 1 -^,2 C ;;+ 1 )-2fcD c (4 /!C ;;+ 1 ,< +1 )+fc((, /1 +A)< +1 , c ;; +1 ) 

+fcb,K,^ +1 ^ +1 ) = fc(4,< +1 ). 

We have (c£ +1 - c£,2c£ +1 ) = |c^ +1 | 2 - |c£| 2 + | C ;; +1 - c£| 2 . Thanks to propositionsgjandgr 

~2k(A h cl + \cl +i )=2k\\c" h +l \\l 1 b h (u" h ,4 +1 ,c" h +1 )^0. 
Using the Cauchy-Schwarz and Young inequalities we write 

Finally thanks to ( 12.1 U and ( 13.41 ) we have 5/, ^ 0. Thus we obtain 

|c;; +1 | 2 -| c ;;i 2 +2fcA.|| c r 1 ||^+4l c ;; +1 l 2 <"- 
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Let ra £ { 1 , . . . , N} . Summing up the latter relation from n = to m — 1 we get 

m m 

'2kD c ^\\c" h \\U\c° h \ 2 + C^k^C, 



-A I -t^a.^, 

77=1 /7=1 



thanks to proposition ^. II With a similar work on equation ( 15.3b . we get (16.21 1. We now prove (16.31 1. Let 
n = m— 1 £ {0, . . . ,N — 1}. Multiplying equation ( 15.4| > by —p1 +l and using proposition ^. 11 we get 

(^ +1 v^ +1 , v^ 1 ) = (C 1 , v /l7 ,;; +1 ) + (n^pf 1 ). (6.4) 

The left-hand side term satisfies (k"/" +1 V<,/^ +1 , V/,/?)J +1 ) ^ 1 V/,/^ +1 1 2 . We now consider the right- 
hand side. Using 15.11 . the Cauchy-Schwarz and Young inequalities we write 

Also, the stability of fl/>nr for the L 2 -norm, proposition ^. H and the Young inequality lead to 

\(n Pr s h , P t +l )\ < \s h M +1 \ ^c\ P 'i +l \ <c\v hP i +1 \ ^^L\v hP i +l \ 2 +c. 

Thus we deduce from d6T4j» that | V,,^ +1 | 2 = | V,,p™| 2 ^ C. Then d53]l imply 

K\ = K +1 i < iCi + kr'vrf 1 ! < iif|i L . (0ir;L 2 ) + i| f c|| wl . ((0 , 1)x(0H) |v /!/ ,»+ i i ^c. 

Estimate ( 16.3b is proven. ■ 
7. Convergence analysis 

Let e = max(/i,A:). In this section we study the behavior of the scheme (I5.2b -( l5.5b as e — * 0. We first 
define the applications c £ :R ^ P , c £ :R ^ P Q , 9 £ :R ^ P Q , p £ :R^ Pf c , s £ : R -> P , 4 : K -> -Po 
and u E : K -> RT , ff : M -> P by setting for all n £ {0,...,N— 1} and t £ [r n ,f„+i] 

c e (0 = CJ+ 1 , c e (0 = c\ + I (r -f„) (cf 1 -cj), e (O = 
Pe(t)=pt + \ s e (t)=s h , s e (t)=s c h , u e (t)=u n h , i e {t)=^\ 

and for all t <£ (0, T) 

c E (t) = c £ (t) = 9 £ (t) = p £ {t) = s £ {t) = s L £ {t) = 0, u £ (f) = f £ (f) = 0. 

We recall that the Fourier transform / of a function / £ L 1 (R) is defined for any T £ R by 

/(t)= / e- 2i ™f(t)dt. (7.1) 



We begin with the following estimate. 

1 

4' 

"|Tr(|c e (T)| 2 + |0 e (T)| 2 )rfT<C. 



Proposition 7. 1 Let < y < I. There exists C > such that for all e > 
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PROOF. Since equations (15.2b and ( 15.3b are similar we only prove the estimate on c £ . We first define 
g e : K — > Pq n Lg as the solution of 

A ft g e = D c A h c e + s c e - (s e + A ) c £ - b h (u £ , c £ ) . 

Multiplying this equation by — g £ we obtain 

- {Ah8e,ge) = -D c (A h c e ,g s ) - (s c £ - (s £ + X)c £ ,g £ ) +b h (u £: c £ ,g £ ). (7.2) 

Proposition |4.5| allows us to write 

-(A h g e ,g e ) = \\g £ \\ 2 h , -{A h c £l g £ ) ^ \\c £ \\ h \\g £ \\ h . 

Thanks to the Cauchy-Schwarz inequality, d5.11 l and proposition |3.2| we have 

|(4 - (s £ + l)c £ ,g £ ) | < CQs c \ + \s\+X) \g £ \ C\\g e \\ h . 

According to proposition l4.3l then proposition ^. H and (16.31 . we have 

\bh(u e ,c e ,g e ) | s$ C||c E ||o= \\g e \\h |divu £ | +C\\c £ \\ h \\g £ \\ h \u £ \ 
^C\\g £ \\ h \dwu £ \+C\\c £ \\,,\\g £ \\ h . 

Let us plug these estimates into ( 17.21 ) and integrate from to T . We get 

\\g £ \\ h dt <C / r |divu £ |A + C f \\c s \\ h dt^C, 



because of proposition 15.11 and ( 16.2b . Definition (17. lb then leads to 

VreM, |||e(T)|U<C. (7.3) 

We now use this estimate to prove ( 17.1b . Equation ( 15.2b reads 

^c £ = A h g £ + (c° h So - 5 7 ) 

where 5a and 5t are Dirac distributions respectively localized in and T. Let Tel. Applying the 
Fourier transform to the latter equation we obtain 

-2mt5 e (r) = A h g e {z) + (c° h - c»e- 2 '™ T ). 

Let us take the scalar product of this relation with /sign(r)c £ (T). Applying propositions 13.21 and 14.51 
leads to 

2tt|t| |? £ (t)| 2 < C (||&(t)|| a + \c%\ + |c£|) \\5 e (r)\\ h . 
We assume that T ^ and multiply this estimate by |t| 2 ' 1 '~~ . Using proposition ^. 1 l and ( 17.31 ) we get 

ItP^wi^cItI^-'ii?,^)!!,. 

Using the Young inequality and integrating over {t G K; |t| > 1}, we obtain 

\T\ 2 r\5 £ (T)\ 2 dT^ I ^^dX + cj \%{x)f h dx. 

T >1 J t >l J t >l 
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For |t| ^ 1, we have |T| 2r |? £ (t)| 2 < |? £ (t)| 2 ^ C||? £ (t)|| / 2 according to propositionEOl Thus 

r\ 2 r\%(x)\ 2 dr<C f \\%(T)\\ 2 h dx. 



|t|<1 -'\t\<\ 
By combining the bounds for |t| > 1 and |t| ^ 1 we get 

/ \T\ 2 ?\5 £ (T)\ 2 dT^ [ \T\'y- 2 dX + C ( fc £ {T)\\ldX. 

Jm. J\r\>\ Jm 

Since 4y — 2 < —1, we have Jj^i |t| 4 ^ — 2 c/t ^ C. Thanks to the Parseval theorem and (16.21 ) 

/ ||? e (T)||fcT< / |^||^<C(*||c»||2 + *£||cg||lUC, 

Jm. Jm V y 

because = ||n/> co||/, < C||co|| 1 (see fl3tl p. 776). Hence the result. ■ 
We can now prove the following convergence result. 

Proposition 7.2 There exists a subsequence of (c e ,8 £ ,p £ , u £ ) £ >o, not relabeled for convenience, such 
that the following convergences hold for e — > 

c £ ->c inL 2 (0,7/;L 2 ), 8 £ -> in L 2 (0, 7/;L 2 ), (7.4) 

p £ -^p weakly inL^O.Jji? 1 ), u £ u weakly in L 2 (0, 7/;L 2 ). (7.5) 

The limits (c,0,p,u) satisfy the following properties. We have c G L 2 (QJ\H X ), 9 € L 2 (0,T;H l ), 
p G L°°(0,T;H l ) and u G L°°(0,r;L 2 ). We also have < c(x,t) < 1 and 0" sC 9(x,t) < 0+ a.e. in 
i2 x [0,T]. For all G ^"(fl X (-1,7")), c and satisfy 

jf ( (c, d, ) + D c ( Vc, V0 ) - c (u • V0 ) - (s c - (s + A ) c) ) dt = (c , ( • , 0) ) , (7.6) 
0,d,0) +D e (V0, V0) - (u- V0) + ( Jfl +s(6 - 9,)) <j>) dt = (0o,0(-, 0)). (7.7) 



Lastly we have 

U = f-Jc(e,0)Vp inL°°(0,7/;L 2 ), divu = s inL 2 . (7.8) 

Proof. In what follows, the convergence results hold for extracted subsequences. They are not rela- 
beled for convenience. We begin by proving ( 17.4b . According to proposition ^. II the sequence (c £ ) £> o 
is uniformly bounded in L°°(0, T;L 2 ). Thus there exists c G L°°(0, T;L 2 ) such that 

c £ — * c weakly in L 2 (0, T;L 2 ). 

Using the Fourier transform, we prove that this convergence is strong. Let d £ = c £ —c and M > 0. We 
use the following splitting 

\d £ (T)\ 2 dr= [ \d £ (T)\ 2 dx+ [ \d e {x)\ 2 dx = I e A +jf. (7.9) 

J\t\>M J\t\^M 
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Since \d £ (%)\ 2 ^ 2|? £ (t)| 2 + 2|?(t)| 2 we have 

/ £ M ^2 / \? £ {T)\ 2 dT + 2 [ \c(T)\ 2 dT. 

J\x\>M J\r\>M 



Hence 



Using proposition l7.1l we write 

/ , \^(t)\ 2 dx < / \r\ 2 7\c e (x)\ 2 dx < 

J|t|>M m z 'J| t |>m M z T 

M-i J\r\>M 

This implies that for all £ > 0, — > when M — > °o. We now consider Let T € [— M,M]. Since 
CeW € for all f G K, and c £ c weakly in L 2 (Q,T;L 2 ), we deduce from {7J) that c £ (t) G Pq and 
c £ (t) — 1 c(t) weakly in L 2 . Extanding c £ (t) by outside i2, one checks ( 11311 p. 81 1) that 

V77GK 2 , |c e (T)(- + ij)-c e (T)|<C||c e (T)|U|Ti|(|T/|+A). (7.9) 

Then, using estimate ( 16.2b . we deduce from fol] (p. 834) that c £ (t) — > c(t) strongly in L 2 . Thus <i £ (t) = 
c £ (t) — c(t) — > in L 2 , so that 7^ — > when e — > 0. Now, let us report the limits for Iff and Jf^ into 
( 17.91 . Using the Parseval identity we get 

2 / j |2 j, / i „|2 



£ / £ (T)| Z flfT= / 141^?= / \c e -c\ Z dt->0. 

JR JR 

Thus we have proven that c £ — > c in L 2 (0, T\L 2 ). A similar work proves that £ — > (9 in L 2 (0,T;L 2 ) 
with G L°°(0,r;L 2 ). Hence {73) is proven. Moreover using proposition 16 . 1 1 we obtain ^ c(x,t) ^ 1 
and 0" < 0(x,f) sj 0+ a.e. in Q x [0,7/]. Lastly, using (J62) and (fL9l . we get as in 0] (p.811) that 
ceL 2 (0,T;H l ) and G L 2 {0,T;H l ). 

Let us now consider the sequences (p £ ) £ >o and (u £ ) £> o- According to ( 13.31 ) and ( 16.3b the sequence 
{TIpcp £ ) £> Q is bounded in L°°(0,T;H ). It implies that there exists p £ L°°(0,T;H l ) such that Upcp e — s - 
/? weakly in L 2 (0 7 T;H l ). Using proposition 13.31 we get p £ — weakly in L 2 (0,T;H l ). Moreover, 
according to ( 16.3b . the sequence (u £ ) £> o is bounded in L°°(0,r;L 2 ) . Thus we have u £ — ^ u weakly 
in L 2 (0,r;L 2 ) with u G L°°(0,r;L 2 ). We check the properties of u. Using a Taylor expansion, the 
Cauchy-Schwarz inequality, and a density argument, we have 

\\k(c, 0) - k(c s , £ )||l2 (ot . iL 2 ) II k||wi--((o,i)x(o,-)) (\\ c ^ c ^\\l 2 (o,t;l 2 ) + W° - ^WlHo.t-l 2 ))- 

Thus, using the strong convergence of the sequences (c £ ) £ >o and (0 £ ) £ >o, we have fc(c £ , £ ) — * k(c, 0) 
in L 2 (0,7/;L 2 ). Since Vi,p £ — Vp weakly in L 2 (0, T;L 2 ), we deduce from this 

K(c e ,e e )V hPe ^K(c,6)Vp weakly in L 2 (0,r;L 2 ). (7.10) 

Now let v G L 2 (0, T; (^ °°) 2 ). According to (J53) we have 

(u £ ,n RTo v) = (f £ - k(c £ , e £ )v hP£ ,\). 

Using proposition l3.3l one checks easily that (f £ ,rjRT v) — > (f,v) and (u £ ,nRT v) — > (u,v) in L l (0,T). 
Using moreover convergence ( 17.101 ) and a density argument, we deduce from this that u = f — )f(c, ) Vp. 
And since divu £ = s £ by proposition l5.ll we also have divu = s. 
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We finally prove that c satisfies ( 17.6b . For all t G (0, T) equation ( 15.2l i reads 

c/ 



r// 



c £ - D t . A ft e e + b /; (Ug , c £ ) = s c e - (s e + X ) c £ . 



Let y G < ^o°°(X2 x (—1,7')) and l///, = TIp^. Multiplying the latter equation by iff/, and integrating over 
[0, T], we obtain 

(^ e '^' ! ) dt ~ Dc Io ( A i> C£ 'V^ dt + J Q b h(u e ,c e ,Yh)dt = (s L £ -(s £ + X)c £ ,\t/ h )dt. (7.11) 



We now pass to the limit e — * in this equation. We begin with the term J b/,(u e ,c e , Yh)dt. We use 
the splitting b(u,c, y) - b/,(u e ,c e , y/j,) = Af +A| +A3 with 



Af = b(u,c,i/f)-b(u e ,c,v), A| = b(u e ,c,v^)- /_ div(cu e ) y h dx , 
Af= / div(cue)v^rfx-b A (ue,c e ,VA)- 
According to definition (I3.11l l 

Af = b(u,c,y) — b(u e ,c,y) = - / c(u — u e ) ■ Vydx. 



in 

We know that cVi/AGL 2 (0,r;L 2 ). Since u £ u in L 2 (0, T;l?) we get $ A\dt — > 0. We now consider 
Af. We have 



Using the Cauchy-Schwarz inequality we get 



Af = / (W— Vft)div(cu e )dx = / (v/- y//,)(u e • Vc + cdivu e )dx 



jf |A||£fc < || V— V3klU-(Ox(0,r)) (H u ell L 2 (0ir . L 2 ) + Hdivuell^^^p My^^iy 

Using a Taylor expansion, one checks that || y- Yh\\v°(nx(0jr)) < /1 ll v Vlk"(flx(0,r))- Thus JoA%dt 
0. Finally we estimate A3. For all triangles K G an d 1 € ^ sharing an edge a, we set c^,l = c# if 
Ue • riK.o ^ and cjc.l = cl otherwise. Using the divergence formula, we deduce from definition ( 13.121 ) 



A 3 = E E Wk {c-c K , La ){\x £ -a K , a )da 
= E (YK a -WL a ) I {c-c Ka<La )(u e -n Kat a)da. 



Using definition ( 13.51 ) this reads 

A3 = E (fc-fj [ {npncC-C Ka ,L a ) (He- n Ka ,<T) dO 



E (Vk - VO \a\ ((npncc)(x a )-CK a ^ (u £ -n Ka , a ) ( 
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Using a Taylor expansion, one checks that \ Yk„ — \ ^ h II^ 7 V / 1Il°°(I2x(o,7'))- Moreover \a\ ^ h. Thus, 
using the Cauchy-Schwarz inequality, we have 

|Af| < Ch 2 £ |u e (x CT )| (n PT c){x a )-c K ^ La 



^ ch 2 ( £ |u £ ( Xo -)| 2 ) 1/2 ( £ 



,2\l/2 



Using the assumption on the mesh, one checks that 1^1 ^ Ch 2 for all K £ 3^,- Thus, thanks to a 
quadrature formula, we have 



K\ 



* 3 



1/2 



Er E |(n Prc )(x a )- c ^| 



1/2 



< C|u e | \TIpncC-Ce\- 

We write np<cc — c £ ~ (Ilpncc — c) + (c — c £ ) and we use proposition l3.3l . We obtain with ( 16.3b 

\Al\dt^C\\u £ \\ L „^ rL 2 ) (h\\c\\ L 2 {0T . Hl] + \\c- c e \\ L 2 (0tT . M2) ). 

since c £ — > c in L 2 (Q,T;L 2 ) when e = max(/i,fc) — > 0, we conclude that / A3 off — > 0. Gathering the 
limits for Aj, A%, A3, we obtain 



b h (u E , c £i n Po \j/)dt — > y b(u,c, 

We now consider the other terms in ( 17.1 It . Proposition |4.6| leads to 

(A h c e , n Po W) = {c e , A h (n Po y/)) = (c e , A A (n Po y) - A y) + (c e , A y) . 
According to proposition |4.7| 

(c £ ,A h (n Po y)-Ay)\ < ||c e || ft \\A h (n^y) - ^^ll-M < ch IMI* Wh- 
We then apply the Cauchy-Schwarz inequality and use ( 16.2b . We obtain 



(7.12) 



£ J (c £ ,A h (n PoV ) -A V )\dt^ch (J* \\c e \\ldt) l/2 ^ch(k£\ 1 4 1 1 1) ' /2 C h 



Moreover, since c £ — > c in L 2 (0, T;L 2 ), we have Jq (c £ ,Ay)dt — > [J (c,Ay)dt. Thus we deduce from 
(I7T21 



We are left with two terms. First, using Taylor expansions, one checks that 

y h ^y, d t y h ^d t y mL 2 {Q. x (-1,7)), l^(-,0) v(-,0) inL 2 . (7.13) 
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We know that c £ -> c in L 2 (0, T;L 2 ). Thus 

/ {s c -(s + X)c e ,y h )dt^> [ (s c -(s + X)c,Y)dt. 
Jo Jo 



Finally, integrating by parts the first term of ( 17.1 j} , we get 

\dt~ 



\ f T 
-c e , Yhj dt = (c e , y/h)t=T - {c e , Wh)t=o - J o {ce,dtYh)dt. 



Since y/ e x (— l,,T)) we have (c E ,Yh)t=T = 0. Using proposition 13.31 one checks that c° 

Ilp Q C() — > cq in L 2 ; using moreover ( 17.131 ) we get 



(c £ ,Vh)r=o = (cl,Yh(-,0)) = (n Po c ,r/,(-,0)) -> (c ,r(-,0)). 

2 (o,: 



For the last term, one easily checks that \\c e — c £ \\ L 2/ 0T . L 2\ — ► 0. Thus, since c £ — > c in L 2 (0,T;L 2 ), 



we also have c e — > c in L 2 (0,T;L 2 ). Using moreover ( 17. 1 3b we get J (c e ,d t \j/i,)dt — > J (c,d t \j/)dt. 
Therefore 

/ iTt d£,Vh ) dt ~( c °'^''°)) (c,d,W)dt. 
By gathering the limits we have obtained in ( 17. lit we get ( 17.6b . The relation ( 17.7b for is proven in a 
similar way. ■ 

8. Error estimates 

We have proven in section |7]that the problem (12.5b — ( 12.9b has a weak solution (c, 9,p,u). From now on, 
we assume the following regularity for this solution: 

c,0 e^{Q,T;H 2 ), c„d t eL 2 (0,T;H l+r )n^(0,T;L 2 ), 

c„,e„eL 2 (0,r;i 2 ), p£tf(0,T;H 2 ), n6^(0,T;H 1+s ), 

with r > and s > 0. We also assume thatf e ^(O^H 1 ). 
8.1 Definitions 

For all f7z € {0, . . . ,Af}, we define the following errors 

eZ c = c(t m )-c';;, ei 9 =c{t m )-ez, 
e% p =p(t m )-p%, e£ u =u(* m )-<. 

We have the following splittings 

in _ in , in in _ in , „m 

e /?,c — fc A, c + '//i.cj e A ; e — fc «,e + w/i.ei 

e h,p — fc /!,p "T f //!,/IJ e /l,U — fc /l,u+ '//l,0 



with the discrete errors 



e h,c - n Po c(t m ) - c™, e"' e - n Po o (t m ) - e™, 
e)l P = n Pr p(t m ) - pz e »; u = n RTo u(/ m ) - <, 
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and the interpolation errors 

< c = c (f>n) - n Po c{t m ), 77;;; = e(t m ) - n Po e(t m ), 
nZ P = p(tm) - npnc P (t m ), 77/;; u = u(t m ) - n RTo u(/ m ). 

The interpolation errors are estimated as follows. We write |tj!",| ^ \c(t m ) — np Q c(t„,)\ + \Ilp Q c(t m ) — 
np c(t m )\ and the same for r]™ d . Using proposition l3.3l and ( 13.7b we obtain 

\nt\< ch M^)h< ch \\ c \\ L -m^)> \riZe\<ch\\e\\ Ir(0tTiBl) . (s.i) 

According to proposition |3.3| we also have 

K P \ + |V/«C P | ^ Ch \\p{t m )\\ 2 ^ Ch |b||^ (0ir;ff 2 )5 (8.2) 
|77, m u |<O||u(? m )||i<C/ J ||u|| L » (0!r . Hl) . (8.3) 

We now have to estimate the discrete errors. 

Proposition 8. 1 For all n G {0, . . . ,N — 1} and \j/ h G P" c we have 



k 



e - D c A h e^ + b h (e» n Pa c(t n+l )) + b h (ule" h + l ) + (s" h + A ) = Cft 1 + C£+\ (8.4) 



£hM ~ £hd -Deleft? +b h ( £ '^n Po e(t n+l )) +b h (u" h ,eZy)+s h e$ = + (8.5) 

(K{cl + \d^)V h e^\V hWh ) =-((^ 1 e^ 1 + ^ 1 <+ 1 )Vp(Wi),V ft% ) - (Pr\V ftW ,), (8.6) 

e^-i^Kt^ + ^l^ (8.7) 

For all m G {0, . . . ,#}, the consistency errors q" p q" 2 , ©£\, ®j™ 2 , P£ ! and UJ are defined in ( fOt , 
( 18.101 ), ( 18.141 ). ( 18.15b and the terms JC^j and fc™ 2 are given by ( 18.131 ) below. 

Proof. Let n G {0, . . . ,N- 1}. Equation ( TO] ) for f = ?„ + i reads 

^c(f„+i)-Z) c 4c(f n+ i)+b(u(/„ + i),c(f n+ i)) = s c -(s + k)c(t n+ i). 
We introduce the time discretization by setting 

Rn+i = ( c{t n+] )-c{t n ) _ Ct{tn+i) ^ + b( u ( fn )-u( f „ + i),c(Wi))- 

We get 

c{t n+1 )-c{t n ) _ DcAc{fn+i) + b(u( / „), c fa + i)) = *c - (* + A)c(Wi) 
We apply flp to this equation. By subtracting the result from ( 15.21 ) we get 

+iT n) b(u(f n ),c(f n+1 )) -b/^u^cf 1 ) + i7 / . (( i + A)c(r„ +1 )) - (^ + A) C ;; +1 -n Po ^+i. (8.8) 
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We now introduce the discrete errors as follows. Since c(t n+ \) — c(t n ) = J^" +1 c t (s)ds one checks that 
„ ( c(t„ +l )-c(t n ) \ c';+ l -c n h l['"^ (r[r( , ff ^MU.x 1 ^! r» \ 

We also have 

n Po Ac(t„ +l ) -A h c n h +1 = n Po Ac(t„ +1 ) - A h (n Po c(t„ +1 )) + A h e^+\ 
Using the linearity of b/,, one easily checks that 

n Po b ( U (r„) , c(t n+l ))-b h (»ic>;+ 1 ) = b h (ui ) +b h (efc n Po c(tn + i)) 

+ np g b (u(t„) , c(t n+ 1 )) - b h (n RTo u(t„),n Po c(t n+ i)). 

Lastly 

n Pa ((s + l)c(t n+i )) - {s h + X)cl +l =77p ((, + A)7 7/ :+ 1 ) + {s h + X)el+\ 

Using these relations in ( I8.8I 1 we get (18.41 1. For any m G { 1 , . . . , A^}, the consistency errors G Po and 
C™ 2 G Po are given by 

Cft = n f0 ( c(fm) ~ fc c(f '"- 1 } - c f (r m ) + b(u(f m _ 1 ) - ufa) , c(r m )) ) 

1 r'>" ~ 

+ n Po {{s + X)r]l c )~- (n Po c t (s)-n Po c t (s))ds, (8.9) 
k Jt m -\ 

C% 2 = D c (n PQ Ac{t m ) - A h (n fb c(t m ))^j - (ll^b(u(/ m _i),c(/ m )) -b/,(i7 RTo u(f m _i),n / . c(f m ))) . 

A similar proof leads to ( 18.5b where the consistence errors @™j G Po and &]" 2 G Po are defined for any 
m G {1,...,A^} by 

0™i = n Po ( e(fm) -^( f "-') _ 9t (t m ) +b(n(f w -i) - n(f„), flfa))) 
1 /■'"> ~ 

+ n Po (sT]Z e )-- (n Po d t (s)-n Po d t (s))ds, (8.10) 
© A m 2 = D e (n Po 40fa)-4 A (n Po 0fa))) - (n Po b(u(t m ^),e(t m )) -b h (n RTo u{t m ^),n Po e(t m ))y 

We now consider the problem associated with the pressure. Let n G {0, . . . ,N — 1} and l//), G P" c . Mul- 
tiplying equation ( 12. 7t written for f = f„ + i by y/h ar >d integrating by parts, we get 

(jc(c(/ B+ i),0(/„ + i))Vp(?„+i),V A v^) = (%+i)>VftV^) + (^V^)- (8-H) 
On the other hand, using ( 15.4b and proposition l4.ll we have 

Since V,,^/, g P , one checks that (f£ + \ V A y/j,) = (JT Po % +1 ), V A y/i) = (f(f„+i), V/^/,). According to 
(13.5b we also have (LZ/vks/,, y// 7 ) = (s/,, y//,). Thus 

(k(c£ +1 ,0 v^ +1 ,v ftWl ) = (ffwO.Vtfti-is/nft). 
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Substracting ( 18.1 U from the latter relation, we obtain 

(K(c(t n+l ),0(t n+l ))Vp(t n+1 )-K(4 + \eZ +l )V h pl + \V h y h ) = -( S - Sh ,y h ). (8.12) 
We split the left-hand side as 

K{ci + \e^){yp(t n+l )-y hP i +l ) + {K{^ 

Using a Taylor expansion, one can check that 

K(c(t n+l ),e(t n+1 )) K(cl + \e» +l ) = (e^ 1 + O *-/u 1 + Ke + ^)Kf- 

We have set for any m 6 {0, . . . ,N} and s G [0, 1] 



We also have 



<%{ S ) = ci + (c(t m ) -c%)s, e?(s) = e;: + (e(t m ) - e?)s, 

Ka = f K x (c%(s),8?(s))ds, k™ 2 = f K y (c>»(s),e™(s))ds. 

Jo ' Jo 



(8.13) 



Plugging these relations into (I8.12l i we get ( 18.6b . For all m e {0, . . . ,N} we have 

K = + <2<e) Vp(^) + *(« m ) (8.14) 

We end with the equation associated with u. Let n S {0, . . . ,iV— 1}. Applying the operator /Trt,, to 
( I2.7I > for t = t„ + i we obtain 

rt R T u(f„ + i) = i7 R T f(f n +l) -iT R T o (K-(c(f n+ i),0(f„ + i))V/7(f n+1 )). 

Let us substract this equation from ( 15.51 ). Since f? + = rip f(f, 7+ i) we get 

n RTo u(f„ +1 )-uf 1 = 77 RTo (f(f„ +1 )-n Po f(f„ +1 )) 

- n RTo (k (c(t n+l ),e(t n+l )) v P (t„ +l )- < +1 v /)P ;; +1 ) . 

One easily checks that 

iT R T u(f„+ 1 ) - u;; +1 = 77 RTo (u(f„ + , ) - iT RTo u(f n+ , ) ) + e£+ 1 
Thus we get (18.7b . For all m € {0, . . . ,N}, we have 

u;; j = n RTo ((f (t m ) - n Po t (t m )) - n »- p" h ') . (8.15) 

This ends the proof of proposition l8.ll ■ 
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8.2 Error estimates 

We first estimate the consistency errors. 

Proposition 8.2 For all m e {1,. . . ,N} the consistency errors satisfy 

m m 

kj:\cL\ 2 +kj:\®;;,i\ 2 ^c(h 2 +k 2 ), 



n=\ 



n=l 



kl,\\Ci 2 \\\ h +kY,\\®U-i^Ch 2 , 



n=\ 

\n\+\m\^ch. 



PROOF. Let n g { 1 , . . . , A^}. Since the operator Tlp Q is stable for the L 2 -norm we have 
\n Po R"\ < \R"\ < C{tn) ~f"~ x) -c t {t„) + b(u(t n ^)-n(t n ),c(tn)) 



Using a Taylor expansion and the Cauchy-Schwarz inequality, we get 

c(t n ) - c(t„-i) , , .1 f'»-' 



-c t (t n ) ^jf" ' \t„-i~s\\c [t (s)\ds^y^^" ' \c n (s)\ 2 ds) 



1/2 



On the other hand, since Vc(t n )\gQ = 0, we deduce from ( 13.101 by integrating by parts 

b(u(f„_ i ) - u(t„) , c(t„) ) = (u(f„_i ) - u(r„ ) ) • Vc(f„ ) . 
Using a Taylor expansion and the Cauchy-Schwarz inequality, we get 

\u t {s)\ 2 ds 



(u(/„_i)-ufe))-Vc(/ n )|< Vfc|klL~( ,r;tfi) 



Thus 



|ITV«»| ^^{j^ \c,{s)\ 2 ds) 1 +Vk\\c\\ L ~ mHl) [J^ \u,(s)\ 2 ds 
Thanks to the stability of TIp Q for the L 2 -norm and to d8.lt we have 

°°(0,T;L 2 ) °°(0,T;H 1 )- 

The Cauchy-Schwarz inequality and (I3.71 i allow us to write 

|n it c f (5)-n it c r (i)|d5<c/zV^(r m HQWIIi+r^ 1 



1/2 



By plugging these estimates into definition ( 18. 9t we get 

k\Cl x \ 2 ^ k 2 \\c\\ 2 L ^ Hl J'"^ \xx,{s)\ 2 ds + k 2 J''^ \c tt {s)\ 2 ds 

2 

\l°°(0,T;H 1 )- 



Ch 2 / " \\c t {s)\\\ +r ds + Ckh 2 \\c\\ 2 



(8.16) 

(8.17) 
(8.18) 



<B-1 
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Summing up the latter relation for n = l to m € {1,... ,N} and using a similar work on ©"j we get 
( 18.16l >. Now let 77 G {1, . . . ,N}. Using propositions l4.7l and l4.4l we have 

1 1 n Po Ac(t n ) - A h (n Pa c(t n )) 1 1 _ i j, ^ Ch 1 1 c\ | £ - (0) r^a) 

and 

\\n Po b{u(t n ^),c(t n ))-b h {n RTo n{t„^i),np c(t„))\\^ij, ^Ch\\c\\ L ~ {0J . H[) \\u\\ L ~, {0T . ll i +s) . 

Plugging these estimates into definition (|8.9l l and summing up from n — 1 to m, we obtain 

m 

k H ll^,2ll-l,A ^ Ch (W C Wl°°(Q.T:H 1 ) H U Hl-(0,7';H 1 + s ) + W C Wl°°(0,T;H 2 )) ■ 

A similar work on @ ; ™ 2 then leads to (18. 17b . We finally prove (18.18b . Let we {1, . . . ,N}. On the one 
hand, we have by ( 18.14b 

W\ < (IKill-TO + ll^ll-Kel) |Vp(f m )l + WKic^e^Uv^j. 

Using estimates d8. lb — d8.31> we get 

l P /"l «S Ch ||K-ilwl-((OJ)x(0,oo))(il C lli"(OT;//') + II \\l~(0,T;H^ ) + WpWu-frT-Jfi))- 

On the other hand definition (18.15b leads to 

\v'i;\ ^ |n RTo (f(^)-JTp f(^))| + |<l+|iT|. 

Using the stability of I1rt for the L 2 -norm and proposition 13. 3 1 we have 

/TR To (f(/m)-n Po f(/ m ))| < \f(t m )-n Po t(t m )\^ch\\t\\ L ^ OJ . Ml) . 

Using moreover ( 18.3b we obtain 

|U™| < Ch\\K\\ w Lo.( {01 ) x ( ^(\\c\\ L °°( 0T . H i- ) + ||0||L-(O,r;//>)) 
+ C{\\p\\ L ~(oj-tf-) + \\f\\[J"(O.T;H l ) + IMIi^O.r;!! 1 ))- 

We have proven (18.18b . ■ 
Using the former proposition we are now able to estimate the discrete errors. 
Proposition 8.3 There exists some real ko > such that for any k < ko and m e {1, . . . ,N} 

m 

Kc\ 2 + Ke\ 2 + k^ (\K c .\\l + \K e \\l) ^C(h 2 + k 2 ), (8.19) 

n=l 

\V,K P \ + Ku\^C(h + k). (8.20) 
Proof. Multiplying d8~4l i by 2/ce / "+ 1 , we obtain 

g«+l gf! 

( h ' c k h ' c ,2keif) -2D,fc(4 /Ie ;;;^ e ;;; 1 )+2fcb /? (u^ e ;'.+ 1 , e ;;,+ 1 ) +2U \e^ l \ 2 

= 2k(C£{+CZ+\e«+ 1 )- 2k(s h , \e>^\ 2 )-2kb h (e^n Po c(tn +l ),e" h + [ ). (8.21) 
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Using an algebraic identity we have 



2Jt| £l) = | £ »+t |2 _ |C| 2 + |£ „+1 _ g „ c| 2 



We know by propositions 14. 2l and l4.5l that 

-2k(A h e£\eZ+ 1 )=2k\\e£X, b*^,^ 1 ,^ 1 ) >0. 

We have 

2fc(, A ,| eft "+ 1 | 2 )<2A:|^|| I>s |6^ 1 | 2 <CA:K+ 1 | 2 . 
Using the Young inequality, we also write 



■-I 

2 

We are left with the termb/,(£^ u ,n/> c(f„ + i),e / "^ 1 ). We have £" 2U <E RTo. Using the divergence formula, 
proposition l5. 1 l and (13.5b . one easily checks that dive^ u = 0. Thus we can apply proposition l4.3l to get 



|b /I (e,« l( n / . c(f„ +1 ),<+ 1 )| <ci^j ||n^ c (r„ +1 )|| A H^+'IU- (8-22) 

Let us first bound ||n/> c(f„ + i)||/,. We have 

\\n Po c(t, 1+l )\\ h s$ \\n Po c(t„ +l ) - n Po c(t n+l )\\ h + \\n Po c(t n+1 )\\ h . 

Using an inverse inequality (see proposition 1.2 in I22I0 and ( 13.7b 

~ C ~ 

\\n Po c(t„ +l ) - n Po c{t n+l )\\ h ^ -\n Po c(t„ +l ) - n Po c(t„ +l )\ ^c. 

n 



Moreover, according to 1131] (p. 776), we have \\n Po c{t n+x )\\ h ^ C\\c\\ L «, {QJ . Hl) . Thus \\n Po c(t n+l )\\ h ^ 

C. We now estimate |£jJJ. Using the stability of /Trt for the L 2 -norm and the Cauchy-Schwarz in- 
equality in ( 18.7b we get 



KJ < IMIwi,~((0,i)x(0,~)) (ll V ^ll^ ( o,r;L 2 ) (IC I + KVD + I^C IJ + l U ft I- (8-23) 
E/" p ' as follows. Setting \jf h = e^ 1 



We bound e ; "+' as follows. Setting yff, = e^ +l in (18.6b and using the Cauchy-Schwarz inequality, we get 



(K(cf 1 ,0r 1 )V, e »+ 1 ,V A% "+ 1 ) < l|jf||^((o,i)x(o > -))ll V Pll i -(o I ^ 2 )(l e V 1 l + l<Vl)|V*< 1 | 



+ iprnv^i+k-^iKp 



The left-hand side is such that (K(t^ +1 ) 0^ +1 )V A e^+ 1 ,V A ^+ 1 ) > JQ n/ 1 V/.e^ 1 1 2 . As for the right-hand 
side, we have Sf, = n P(y s and e£t € f" 6 iHLq. Thus, according to propositions [XT] and [33] 

k-jAll^KCAWi-^ijiv^ 1 !. 
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Finally |P); +1 1 sC Ch thanks to ( EE) , Therefore we obtain 

\Vh£Zf\<C{h + \elf\ + \el$\). (8.24) 
Let us plug this estimate into ( 18.23b . Since |UJJ +1 1 < Ch thanks to (18.18b . we get 

I^J^C^ + le^'l + le^ 1 !). (8.25) 
Now, plugging this bound into (18.23b and using the Young inequality, we obtain 
k\b h {el^c{t n+l ),elf)\ < Ck{h + \elf\ + \e$\) \\e"+ l \\ h 

< Dc^\\e^\\l + Ck(h 2 + \£^\ 2 + \e^ l \ 2 ). 
Now we have treated all the terms in ( 18.211 ). This equation implies 

le^+M 2 - |e^.| 2 + Z) t . J t|| e »+ 1 !li < C^(^ 2 + le/"^ 1 ! 2 ^ le^ 1 ! 2 ^ I^/",! 1 ! 2 ^ II^ 1 !!^!./,)- 
Let m £ { 1 , . . . ,N}. Let us sum up the latter estimate from n = to m — 1 . Thanks to (13.3b 

= l^o c <>-c£| = \n Pa co-n Po co\ < Ch\\c\\ L ~ mH2) . 

Using moreover estimates ( 18.161 ) and d8.17| ) we get 



l< c | 2 + J D^Lll<cllNc^(| £ Li 2 +| e ;;, e i 2 )+c(/ ! 2 +^). 

n = 1 11= 1 

Summing this relation with the one obtained by a similar work on < |8.5t we obtain 

m m 

ICI 2 + Ke\ 2 + k£ (D c KJ 2 h +De \\e" Ke \\ 2 h ) ^Ck£ (|< t .| 2 + |< e | 2 ) + C(/ ! 2 + k 2 



Using a discrete Gronwall lemma (see lemma 5.2 in l22ll ) we get (18.19b . Then (18.24b and ( 18.251 ) imply 
(HQ). ■ 

By combining proposition l8.3l with estimates ( I8.1b - (l8.3b . we obtain finally the following result. 

Theorem 8.1 There exists a real ko > such that for all k < ko and me {I,... ,N} 

m 

^^^^^(II^JlI + II^Mll^^C^+fc 2 ), 

n=l 

\V h eZ P \ + \e'H u \^C(h + k). 
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